Demonstration: Reading Data

This notebook explains the data output from the GEODYN Starlette-SLR run.

First, load the data using the reading functions.

[1]:
%load_ext autoreload
%autoreload 2

import pandas as pd
import numpy as np

##################################
# INPUT PARAMETERS:
##################################
sat_file = 'starlette'
arc1 = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/analysis/starlette_analysis/'
SAT_ID = 7501001

##################################
# PATH TO DENSITY MODEL RUN of Choice:
##################################
msis86_model = 'msis86'
path_to_msis86 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+msis86_model+'/'


import sys
sys.path.insert(0, '/data/analysis/util_funcs/py_starlette/')
from a_ReadStarlette import ReadStarlette

AdjustedParams, Trajectory, Density, Resids = ReadStarlette(arc1,
                                                           sat_file,
                                                           grav_id,
                                                           local_path,
                                                           path_to_msis86)
The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/goco05s_msis86/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/goco05s_msis86/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/goco05s_msis86/DENSITY/st030914_2wk.goco05s

 Loading data...

Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded

Adjusted Parameters

  • The Adjusted Parameters are read in from the iieout file (Unit 6).

  • These include the State Vector at the start of the ARC, drag coefficient (CD), the solar radiation coefficient (CR), and the adjusted empirical accelerations (GA).

    • the CD is controlled/output by the DRAG card.

    • the GA is ouput by the ACCEL9 card (the 9 parameter general acceleration)

  • AdjustedParams is a dictionary in which the first dimension stores each geodyn run iteration:

[2]:
AdjustedParams.keys()
[2]:
dict_keys([1, 2, 3, 4, 5])

And the 2nd dimension stores all of the adjusted parameters:

[3]:
AdjustedParams[5].keys()
[3]:
dict_keys(['0CD', '0GA', '0GA_t2', '0XPOS', '0YPOS', '0ZPOS', '0XVEL', '0YVEL', '0ZVEL'])

Each adjusted parameter then has the associated statistics:

APRIORI  VALUE
PREVIOUS VALUE
CURRENT  VALUE
TOTAL DELTA
CURRENT DELTA
APRIORI SIGMA
CURRENT SIGMA
[4]:
AdjustedParams[5]['0XPOS']['CURRENT_VALUE']
[4]:
-1228610.639773413

In the case of CD and GA, there are multiple values for each, and an additional dictionary exists for each value

[5]:
print(AdjustedParams[5]['0GA'].keys(),'\n')

print(np.squeeze(list(AdjustedParams[5]['0CD'].keys())))
dict_keys(['0GA 9P 11', '0GA 9P 12', '0GA 9P 21', '0GA 9P 22'])

[Timestamp('2003-09-14 08:00:00') Timestamp('2003-09-14 16:00:00')
 Timestamp('2003-09-15 00:00:00') Timestamp('2003-09-15 08:00:00')
 Timestamp('2003-09-15 16:00:00') Timestamp('2003-09-16 00:00:00')
 Timestamp('2003-09-16 08:00:00') Timestamp('2003-09-16 16:00:00')
 Timestamp('2003-09-17 00:00:00') Timestamp('2003-09-17 08:00:00')
 Timestamp('2003-09-17 16:00:00') Timestamp('2003-09-18 00:00:00')
 Timestamp('2003-09-18 08:00:00') Timestamp('2003-09-18 16:00:00')
 Timestamp('2003-09-19 00:00:00') Timestamp('2003-09-19 08:00:00')
 Timestamp('2003-09-19 16:00:00') Timestamp('2003-09-20 00:00:00')
 Timestamp('2003-09-20 08:00:00') Timestamp('2003-09-20 16:00:00')
 Timestamp('2003-09-21 00:00:00') Timestamp('2003-09-21 08:00:00')
 Timestamp('2003-09-21 16:00:00') Timestamp('2003-09-22 00:00:00')
 Timestamp('2003-09-22 08:00:00') Timestamp('2003-09-22 16:00:00')
 Timestamp('2003-09-23 00:00:00') Timestamp('2003-09-23 08:00:00')
 Timestamp('2003-09-23 16:00:00') Timestamp('2003-09-24 00:00:00')
 Timestamp('2003-09-24 08:00:00') Timestamp('2003-09-24 16:00:00')
 Timestamp('2003-09-25 00:00:00') Timestamp('2003-09-25 08:00:00')
 Timestamp('2003-09-25 16:00:00') Timestamp('2003-09-26 00:00:00')
 Timestamp('2003-09-26 08:00:00') Timestamp('2003-09-26 16:00:00')
 Timestamp('2003-09-27 00:00:00') Timestamp('2003-09-27 08:00:00')
 Timestamp('2003-09-27 16:00:00') Timestamp('2003-09-28 00:00:00')]
[6]:
AdjustedParams[5]['0GA']['0GA 9P 11']
[6]:
{'APRIORI_VALUE': 0.0,
 'PREVIOUS_VALUE': 5.536419290656419e-11,
 'CURRENT_VALUE': 5.546253541325767e-11,
 'TOTAL_DELTA': -5.54625354133e-11,
 'CURRENT_DELTA': -9.83425066935e-14,
 'APRIORI_SIGMA': 1.0,
 'CURRENT_SIGMA': 1.6930522e-11}
[7]:
AdjustedParams.keys()
[7]:
dict_keys([1, 2, 3, 4, 5])

And the 2nd dimension stores all of the adjusted parameters:

[8]:
AdjustedParams[5].keys()
[8]:
dict_keys(['0CD', '0GA', '0GA_t2', '0XPOS', '0YPOS', '0ZPOS', '0XVEL', '0YVEL', '0ZVEL'])
[9]:
import plotly.graph_objects as go
import numpy as np

Cd_plot = []

labels = list(AdjustedParams[5]['0CD'].keys())

for i in AdjustedParams.keys():

    Cd_plot.append(float(AdjustedParams[i]['0CD'][labels[0]]['CURRENT_VALUE']))

fig = go.Figure(data=[go.Scatter(x=list(AdjustedParams.keys()),
                                   y=Cd_plot,
                               mode='markers',
                                    marker=dict(
                                    size=10,
                                    ),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="Drag Coefficient vs Iteration",
    yaxis_title=r"$C_d \,\,\text{(Unitless)}$",
    xaxis_title="Iteration Number",
    )

fig.show()





[10]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
from plotly.subplots import make_subplots


last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val

labels = list(AdjustedParams[5]['0CD'].keys())


fig = make_subplots(
    rows=1, cols=1,
    subplot_titles=("Time Dependent Drag Coefficient",
                    ))
val_list = []
for i in AdjustedParams[last_iter]['0CD'].keys():
    val_list.append(AdjustedParams[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= 'MSIS 86',
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )


fig.update_yaxes( title=r"$C_D (Unitless)$",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Iteration Number", row=1, col=1)

# fig.update_yaxes( yaxis_title=r"$\frac{m}{s^2}$", xaxis_title="Iteration Number",exponentformat= 'power',row=2, col=1)

# fig.update_layout(
#     title="Along Track Accleration Amplitude vs Iteration",    )
# fig.update_layout(
#     yaxis_title=r"$\frac{m}{s^2}$",
#     xaxis_title="Iteration Number" )


iplot(fig)

Trajectory Data

  • The trajectory data can come from multiple sources. The official output should come from the binary orbit file, but that is not working as of January 2020. Instead, we use the orbit file printout (ORBTVU card) to printout the trajectory data to an ascii file (Unit 10 for keplerian and Unit 08 for cartesian). These are then placed in the KEP_TRAJ and XYZ_TRAJ data directories before the temporary files are purged.

  • The Trajectory structure is a dictionary in which the first dimension represents the different satellites, and the second dimension is a Pandas DataFrame containing the trajectory data from the ORBTVU card in which the first dimension stores each geodyn run iteration:

Indexing the first dimension will show the satellites in the GEODYN run that were saved. They are saved according to their SAT_ID. Since the Starlette-SLR runs only contain the Starlette satellite, this is the only one saved.

[11]:
Trajectory.keys()
[11]:
dict_keys([7501001])

The trajectory data can then be indexed by SAT_ID to show the pandas dataframe containing the data. The stored variables can be found with the .columns method.

[12]:
print(Trajectory[7501001].columns)
Index(['YYMMDD', 'HHMM', 'SECONDS', 'X', 'Y', 'Z', 'XDOT', 'YDOT', 'ZDOT',
       'LAT', 'LONG', 'HEIGHT', 'timeHHMM', 'year', 'month', 'day', 'hours',
       'minutes', 'secs', 'Date'],
      dtype='object')
[13]:
Trajectory[7501001].head()
[13]:
YYMMDD HHMM SECONDS X Y Z XDOT YDOT ZDOT LAT LONG HEIGHT timeHHMM year month day hours minutes secs Date
2 30914 105 55.205719 -3902733.45 5158963.35 -3130061.00 -5675.94599 -1433.45008 4725.70967 -25.954786 118.024110 812276.80 0105 2003 09 14 01 05 55 2003-09-14 01:05:55
3 30914 106 15.005953 -4014288.94 5129496.49 -3035836.70 -5591.75463 -1542.86306 4791.45070 -25.120244 118.880387 811968.94 0106 2003 09 14 01 06 15 2003-09-14 01:06:15
4 30914 106 47.006428 -4190987.23 5077312.08 -2880862.22 -5450.70994 -1718.31837 4893.42231 -23.758977 120.237767 811607.96 0106 2003 09 14 01 06 47 2003-09-14 01:06:47
5 30914 107 13.806875 -4335433.07 5029307.46 -2748618.90 -5327.95616 -1863.82658 4974.68198 -22.607839 121.350810 811436.43 0107 2003 09 14 01 07 13 2003-09-14 01:07:13
6 30914 107 47.207484 -4510756.72 4964053.31 -2580843.47 -5169.22701 -2043.14535 5070.56125 -21.160247 122.709738 811390.75 0107 2003 09 14 01 07 47 2003-09-14 01:07:47
[14]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline


#                     )
x = Trajectory[7501001]['Z'].astype(float)[:]
y = Trajectory[7501001]['X'].astype(float)[:]
z = Trajectory[7501001]['Y'].astype(float)[:]
fig = go.Figure(data=[go.Scatter3d(x=x,
                                   y=y,
                                   z=z,
                               mode='markers',
                                    marker=dict(
                                    size=2,
                                    opacity=0.7),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="3D Cartesian Orbit",
    xaxis_title="X Coordinate",
    yaxis_title="Y Coordinate",
#     zaxis_title="z Coordinate",
    )

# fig.show()
iplot(fig)

Density Data

  • GEODYN does not actually have a built in method for printing out the density. As such, the density has been printd by modifying the DRAG.f90 subroutine in the IIE/MODS directory. The modification consists of printing the density and trajectory data at each timestep of the GEODYN run.

  • The density data is stored as an ascii file named fort.99, which is then stored to the DENSITYdirectory before the temporary files are purged.

  • The density data saves the longitutde, latitude, altitude coordinates, and cartesian coordinates.

[15]:
Density.columns
[15]:
Index(['Elapsed Secs', 'YYMMDD', 'HHMMSS', 'Lat', 'Lon', 'Height (meters)',
       'rho (kg/m**3)', 'delta_rho (kg/m**3/m)', 'X', 'Y', 'Z', 'XDOT', 'YDOT',
       'ZDOT', 'timeHHMMSS', 'year', 'month', 'day', 'hours', 'minutes',
       'secs', 'Date'],
      dtype='object')
[16]:
Density.head()
[16]:
Elapsed Secs YYMMDD HHMMSS Lat Lon Height (meters) rho (kg/m**3) delta_rho (kg/m**3/m) X Y ... YDOT ZDOT timeHHMMSS year month day hours minutes secs Date
0 210.0 30914 434 45.6973 278.5740 1083058.594 1.623031e-15 0.0 182992.974 -5215697.130 ... -1841.620957 -1912.851416 000434 2003.0 9.0 14.0 0.0 4.0 34.0 2003-09-14 00:04:34
1 225.0 30914 449 45.3666 279.6052 1084271.860 1.606673e-15 0.0 284238.712 -5242756.035 ... -1766.174236 -1989.534319 000449 2003.0 9.0 14.0 0.0 4.0 49.0 2003-09-14 00:04:49
2 240.0 30914 504 45.0256 280.6232 1085455.504 1.590705e-15 0.0 385422.960 -5268680.634 ... -1690.382976 -2065.747643 000504 2003.0 9.0 14.0 0.0 5.0 4.0 2003-09-14 00:05:04
3 255.0 30914 519 44.6746 281.6279 1086609.356 1.575121e-15 0.0 486523.867 -5293465.880 ... -1614.263040 -2141.476196 000519 2003.0 9.0 14.0 0.0 5.0 19.0 2003-09-14 00:05:19
4 270.0 30914 534 44.3139 282.6193 1087733.252 1.559916e-15 0.0 587519.621 -5317106.961 ... -1537.830323 -2216.704901 000534 2003.0 9.0 14.0 0.0 5.0 34.0 2003-09-14 00:05:34

5 rows × 22 columns

[17]:


data_nums = 1000


import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density['Date'][:data_nums],
                                 y=Density['rho (kg/m**3)'][:data_nums],
                                    mode='markers',
                                    marker=dict(
                                    size=2,
                                    ),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="Density along Starlette Orbit (for 4 minutes of first day)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)

iplot(fig)








[18]:
data_nums = 10
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density['Date'][::data_nums],
                                 y=Density['rho (kg/m**3)'][::data_nums],
                                    mode='markers',
                                    marker=dict(
                                    size=2,
                                    ),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="Density along Starlette Orbit (Every 10th datapoint)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)

iplot(fig)






Residual Data

  • Similar to the trajectory data, the official output of the residuals should come from the binary residual file (Unit 19), but that is not working as of January 2020. Instead, we read in the residuals that are printed in ascii formats to the iieout file using the OBSVU card.

  • The residuals are saved as a Pandas DataFrame.

[19]:
Resids
[19]:
index YYMMDD HHMM Sec-UTC-R Observation Residual Ratio to sigma Elev1 Elev2 OBS No. ... timeHHMM year month day hours minutes secs microsecs Date Ratio_to_sigma_fixed
0 0 30914 105 55.208831 933.052599576 -0.001484 -0.0015 58.485935 58.485839 1 ... 0105 2003 09 14 01 05 55 208831 2003-09-14 01:05:55.208831 -0.0015
1 1 30914 106 15.009306 1005.251863666 -0.000437 -0.0004 51.260530 51.260437 2 ... 0106 2003 09 14 01 06 15 009306 2003-09-14 01:06:15.009306 -0.0004
2 2 30914 106 47.010250 1145.739853454 -0.002657 -0.0027 41.550447 41.550361 3 ... 0106 2003 09 14 01 06 47 010250 2003-09-14 01:06:47.010250 -0.0027
3 3 30914 107 13.811144 1279.798448021 -0.006932 -0.0069 35.007633 35.007553 4 ... 0107 2003 09 14 01 07 13 811144 2003-09-14 01:07:13.811144 -0.0069
4 4 30914 107 47.212357 1461.073729804 -0.011274 -0.0113 28.408244 28.408171 5 ... 0107 2003 09 14 01 07 47 212357 2003-09-14 01:07:47.212357 -0.0113
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
3999 5688 30927 2306 16.259539 1260.651605954 0.006996 0.0070 53.959755 53.959671 4000 ... 2306 2003 09 27 23 06 16 259539 2003-09-27 23:06:16.259539 0.0070
4000 5689 30927 2306 46.014220 1374.565632664 0.002611 0.0026 46.639968 46.639884 4001 ... 2306 2003 09 27 23 06 46 014220 2003-09-27 23:06:46.014220 0.0026
4001 5690 30927 2307 14.118743 1496.430957455 -0.001994 -0.0020 40.669100 40.669018 4002 ... 2307 2003 09 27 23 07 14 118743 2003-09-27 23:07:14.118743 -0.0020
4002 5691 30927 2307 46.624060 1650.422940823 -0.006786 -0.0068 34.755223 34.755145 4003 ... 2307 2003 09 27 23 07 46 624060 2003-09-27 23:07:46.624060 -0.0068
4003 5692 30927 2308 5.727222 1745.985693706 -0.003580 -0.0036 31.700982 31.700906 4004 ... 2308 2003 09 27 23 08 5. 27222 2003-09-27 23:08:05.027222 -0.0036

4004 rows × 21 columns

[20]:
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=pd.to_datetime(Resids['Date']),
                                 y=Resids['Residual'].values.astype(float),
                                    mode='markers',
                                    marker=dict(
                                    size=4,
                                    ),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="Observation Residuals from Final Iteration",
    yaxis_title="Residuals",
    xaxis_title="Date",
    )

iplot(fig)